load("data/icews95-20.Rdata")
library(dplyr)

df_icews <- mergetable %>%
  dplyr::rowwise()%>%
  dplyr::mutate(source_gov = ifelse(isTRUE(grepl("Government",
                                              source_sectors,
                                              ignore.case =TRUE)==TRUE), 1, 0),
                target_gov = ifelse(isTRUE(grepl("Government",
                                                 target_sectors,
                                                 ignore.case =TRUE)==TRUE), 1, 0))

### filter out gov only; make sure source and target are not the same country
df_icews <- df_icews %>%
      dplyr::filter(source_country != target_country)%>%
      dplyr::filter(source_gov == 1 | target_gov == 1)

## select variables to reduce the size of data file
df_icews <- df_icews %>%
    dplyr::select(-city, -district, -province, -latitude,
                  -longtitude, -story_id, -publisher)

#make year-month variables for aggregation later 
df_icews <- df_icews %>%
      dplyr::mutate(year = format(event_date, format = "%Y"),
                    yearmonth = format(event_date,format = "%Y%m"))


#save for future use
save(df_icews, file = "data/df_icews.RData")
# dim(df_icews)
#[1] 5927685      17
load("data/df_icews.RData")
library(countrycode)

icews <- df_icews %>% 
    dplyr::mutate(source_iso = countrycode(source_country, "country.name", "iso3c"),
                  target_iso = countrycode(target_country, "country.name", "iso3c")) %>% 
    dplyr::filter(!is.na(source_iso)) %>% 
    dplyr::filter(!is.na(target_iso))

# get conflict and cooperation events and their respective score for aggregation
icews_data  <- icews %>% 
    dplyr::mutate(cooperative = ifelse(intensity > 0, 1, 0),
                  conflictual = ifelse(intensity < 0, 1, 0)) %>% 
    dplyr::mutate(cooperative_intensity = ifelse(intensity > 0, intensity, 0),
                  conflictua_intensity = ifelse(intensity < 0, intensity, 0))

icews_data  <- icews_data %>% 
      dplyr::group_by(source_iso, target_iso, year, yearmonth) %>% 
      dplyr::summarise(cooperative = sum(cooperative, na.rm = TRUE),
                       conflictual = sum(conflictual, na.rm = TRUE),
                       cooperative_intensity = sum(cooperative_intensity),
                       conflictua_intensity = sum(conflictua_intensity))
##create the ratio and then take the log
icews_data  <- icews_data %>%  
      dplyr::mutate(relationscores = log((cooperative + 1)/ (conflictual + 1)),
                    relationscores_wt = log((cooperative/cooperative_intensity + 1)/ (conflictual/abs(conflictua_intensity) + 1)),
                    source_ccode = countrycode(source_iso, "iso3c", "cown"),
                    target_ccode = countrycode(target_iso, "iso3c", "cown"))

save(icews, file = "data/icews.RData")
save(icews_data, file = "data/icews_data.RData")

## make it as a matrix

yearmonths <- sort(unique(icews_data$yearmonth))

icews_matrix <- list()
#icews_matrix_wt <- list()

for (i in 1:length(yearmonths)){
  slice <- icews_data[which(yearmonths[i]==icews_data$yearmonth),] 
  slice <- as.data.frame(slice)
  actors <- sort(unique(unlist(slice[c("source_iso", "target_iso")])))
  smat = matrix(0, nrow=length(actors), ncol=length(actors), dimnames=list(actors, actors))
  #smat_wt = matrix(0, nrow=length(actors), ncol=length(actors), dimnames=list(actors, actors))
  # fill in the matrix with matrix indexing on row and column names
  smat[as.matrix(slice[c("source_iso", "target_iso")])] <- slice[["relationscores"]]
  #smat_wt[as.matrix(slice[c("source_iso", "target_iso")])] <- slice[["relationscores_wt"]]
  
  #make it as symetric
  #smat[lower.tri(smat)] = t(smat)[lower.tri(smat)]
  #smat_wt[lower.tri(smat_wt)] = t(smat_wt)[lower.tri(smat_wt)]
  #not define diag
  diag(smat) <- NA
  #diag(smat_wt) <- NA
  
  icews_matrix[[i]]=smat
  #icews_matrix_wt[[i]]=smat_wt
}

names(icews_matrix)=yearmonths
#names(icews_matrix_wt)=yearmonths
save(icews_matrix, file = "data/icews_matrix.RData")
#save(icews_matrix_wt, file = "data/icews_matrix_wt.RData")

library(amen)
## run three AME models
results <- names(icews_matrix)
AME_score <- list()

for (i in 1:length(names(icews_matrix))){
AME_score[[i]] = ame(
  Y = as.matrix(icews_matrix[[i]]), 
  symmetric= FALSE, 
  #nvar=TRUE, 
  rvar=TRUE,
  cvar = TRUE,
  R=1, 
  model='nrm', intercept=TRUE, seed=6886,
  periodicSave = FALSE,
  burn= 5000, nscan= 20000, odens=20,
  plot=FALSE, gof=TRUE, periodicSave=FALSE )
print(i)
}

summary(AME_score)
save(AME_score, file = "data/AME_score.RData")

AME_scorewt = ame_rep(
  Y = icews_matrix, Xdyad = NULL,
  Xrow = NULL, Xcol = NULL,
  symmetric= FALSE, 
  #nvar=TRUE, 
  rvar=TRUE,
  cvar = TRUE,
  R=1, 
  model='nrm', intercept=TRUE, seed=6886,
  burn= 5000, nscan= 20000, odens=20,
  plot=FALSE, gof=TRUE, periodicSave=FALSE )
summary(AME_scorewt)
save(AME_scorewt, file = ".data/AME_scorewt.RData")